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Abstract 

An  implicit  method  for  the  computation  of  unsteady  flows  on  unstructured  grids 
is  presented.  Following  a  finite  difference  approximation  for  the  time  derivative,  the 
resulting  nonlinear  system  of  equations  is  solved  at  each  time  step  by  using  an  agglom¬ 
eration  multigrid  procedure.  The  method  allows  for  arbitrarily  large  time  steps  and  is 
efficient  in  terms  of  computational  effort  and  storage.  Inviscid  and  viscous  unsteady 
flows  are  computed  to  validate  the  procedure.  The  issue  of  the  mass  matrix  which 
arises  with  vertex-centered  finite  volume  schemes  is  addressed.  The  present  formula¬ 
tion  allows  the  mass  matrix  to  be  inverted  indirectly.  A  mesh  point  movement  and 
reconnection  procedure  is  described  that  allows  the  grids  to  evolve  with  the  motion  of 
bodies.  As  an  example  of  flow  over  bodies  in  relative  motion,  flow  over  a  multi-element 
airfoil  system  undergoing  deployment  is  computed. 
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1  Introduction 


Solution  techniques  for  computing  steady  flows  on  unstructured  grids  have  evolved  to  a  high 
degree  of  sophistication.  With  explicit  schemes,  convergence  to  steady  state  is  usually  un¬ 
acceptably  slow,  especially  as  the  problem  sizes  and  complexities  grow.  Therefore,  either 
multigrid  methods  [20,  29]  or  implicit  schemes  [36,  44,  2,  4]  are  required  to  accelerate  the 
convergence.  On  the  other  hand,  solution  techniques  for  dealing  with  unsteady  flows  have 
lagged  behind.  Explicit  schemes,  deemed  to  be  too  slow  for  obtaining  steady  state  solutions, 
may  be  the  schemes  of  choice  for  certain  unsteady  applications,  when  the  time  scales  of 
interest  are  small,  or  more  precisely,  when  they  are  comparable  to  the  spatial  scales.  The 
grids  should  be  clustered  only  in  regions  of  interest;  otherwise,  the  size  of  the  explicit  time 
step  could  become  unnecessarily  small.  However,  when  dealing  with  many  low  frequency 
phenomena  such  as  flutter,  explicit  schemes  lead  to  large  computing  times.  A  time-accurate 
residual  averaging  formulation  [18,  42],  and  temporal  adaptation  techniques  [13],  which  en¬ 
able  different  cells  to  take  a  varying  number  of  local  time  steps  to  get  to  a  particular  time 
level,  can  be  used  to  realize  modest  improvements  in  the  performance  of  explicit  methods, 
but  the  sizes  of  the  time  steps  are  still  controlled  by  the  spatial  resolution.  Therefore,  it  is 
desirable  to  develop  a  fully  implicit  method,  where  the  time  step  is  solely  determined  by  the 
flow  physics  and  is  not  limited  by  the  cell  sizes.  Also,  for  many  practical  viscous  flows,  the 
time  step  restrictions  imposed  by  small  cells  deep  inside  the  boundary  layer  are  excessively 
small.  Since  the  boundary  layer  is  quasi-steady,  implicit  methods  that  allow  for  larger  time 
steps  may  be  more  suitable. 

When  an  implicit  scheme  is  used  to  compute  unsteady  flows,  one  has  to  drive  the  unsteady 
residual  to  zero  (or  at  least  to  truncation  error)  at  each  time  step.  In  the  context  of  factored 
implicit  schemes,  this  is  usually  done  by  employing  inner  iterations  [31,  30,  7,  33].  It  is  the 
role  of  these  inner  iterations  to  eliminate  errors,  if  any,  due  to  factorization  and  linearization, 
and  sometimes  also  errors  arising  from  employing  a  lower  order  approximation  on  the  implicit 
side.  The  number  of  inner  iterations  required  may  be  large  depending  on  the  flow  situation, 
the  size  of  the  time  step  employed  and  the  extent  of  mismatch  of  the  explicit  and  implicit 
operators.  Jameson  [11]  has  advocated  the  use  of  a  full  approximation  storage  multigrid 
procedure  as  a  driver  for  a  fully  implicit  scheme  when  using  structured  grids.  The  significant 
advantage  of  the  approach  when  multigrid  is  used  to  solve  the  nonlinear  problem  is  that 
it  incurs  no  storage  overheads  plaguing  traditional  implicit  schemes  based  on  linearization. 
The  method  is  therefore  particularly  attractive  for  unstructured  grid  computations  in  three 
dimensions.  The  method  allows  the  time  step  to  be  determined  solely  based  on  flow  physics. 
It  has  been  used  to  compute  two-  and  three-dimensional  inviscid  flows  over  airfoils  and  wings 
[11,  1]  using  structured  grids.  Vassberg  [41]  has  applied  this  technique  to  compute  flow 
solutions  over  oscillating  airfoils  using  unstructured  grids  where  a  sequence  of  triangulations 
wa.s  generated  by  removing  points  from  the  fine  grid  triangulation. 

Multigrid  techniques  have  been  successfully  extended  to  deal  with  unstructured  grids 
by  generating  a  sequence  of  non-nested  grids  and  using  piecewise-linear  transfer  operators 
[20,  29].  An  important  development  in  the  use  of  multigrid  techniques  for  unstructured  grids 
is  the  agglomeration  multigrid  algorithm  [15,  37,  43]  for  the  two-  and  three-dimensional 
Euler  equations.  This  technique  has  since  been  extended  to  deal  with  viscous  flows  [14,  22]. 
The  main  advantage  of  the  agglomeration  multigrid  algorithm  is  that  it  only  requires  the 
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fine  mesh  to  be  generated;  the  coarse  grids  are  automatically  generated  by  fusing  fine  grid 
control  volumes  using  an  efficient  graph-based  algorithm,  resulting  in  a  fully  nested  sequence 
of  coarse  grid  levels. 

In  this  work,  the  agglomeration  multigrid  strategy  is  used  to  solve  the  nonlinear  system 
of  equations  at  each  time  step.  The  nested  property  of  the  agglomeration  approach  enables 
a  straight-forward  treatment  of  problems  involving  moving  meshes.  It  is  also  shown  that  the 
mass  matrix  can  be  inverted  indirectly  during  the  multigrid  process.  The  issue  of  the  mass 
matrix  is  critically  examined  in  both  one  and  two  dimensions.  In  order  to  allow  the  grids  to 
conform  to  moving  geometries,  a  technique  is  proposed  and  tested  that  attempts  to  maintain 
the  validity  and  quality  of  the  triangulation.  Inviscid  flow  over  a  pitching  airfoil  and  viscous 
flow  over  an  impulsively  started  cylinder  are  computed  and  the  results  are  compared  with 
experiments  and  other  computations.  Finally,  an  exploratory  computation  is  carried  out 
that  simulates  the  phenomenon  of  flap  deployment. 


2  Governing  equations  and  discretization 

The  equations  governing  compressible  fluid  flow  in  integral  form  for  a  control  volume  V{t) 
with  boundary  S{t)  are  given  by 

4-  [  Wdv  +  <f  [F(W,n,s)  -  G{W,VW,n)]da  =  0,  (1) 

dt  Jv(t)  Js{t) 

IF  =  [p,pV,pef 
F{W,  n,  s)  =  {V  -  s).nW, 

G{W,VW,  n)  =  [0,  t,t.V  -  q.nY, 

where  t  and  q  are  respectively,  the  stress  and  heat  flux  vectors,  p  is  the  density,  V  is  the 
velocity  vector  with  Cartesian  components  Vi,  e  is  the  specific  total  energy,  n  is  the  outward 
unit  normal  vector  of  the  boundary  S{t),  and  s  is  the  velocity  vector  of  the  boundary.  The 
equations  are  augmented  by  the  equation  of  state,  which  for  a  perfect  gas  is 

p  =  (7  -  l)(pe  -  (2) 

In  the  present  scheme,  the  variables  are  stored  at  the  vertices  of  a  mesh  composed  of 
triangles.  The  control  volumes  are  nonoverlapping  polygons  which  surround  the  vertices  of 
the  mesh.  The  contour  integrals  in  Eqn.  (1)  are  replaced  by  discrete  path  integrals  over  the 
faces  of  the  control  volume  which  are  computed  using  the  trapezoidal  rule.  This  technique 
can  be  shown  to  be  equivalent  to  using  a  piecewise-linear  finite-element  discretization  under 
certain  conditions.  For  dissipative  terms,  a  blend  of  Laplacian  and  biharmonic  operators 
is  employed  [20].  The  Laplacian  term  acts  only  in  the  vicinity  of  shocks  and  is  inactive 
elsewhere,  while  the  biharmonic  term  acts  only  in  regions  of  smooth  flow.  Only  the  Laplacian 
dissipative  term  is  used  on  the  coarse  grids. 

After  applying  the  finite  volume  procedure,  the  following  system  of  coupled  differential 
equations  is  obtained: 

4-{VMW)  +  RiW)  =  0  (3) 

dt 
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Here  W  is  the  solution  vector  over  the  whole  field,  -R(W^)  is  the  residual  vector  approximating 
the  second  integral  in  Eqn.  (1),  V  is  the  cell  volume  cissociated  with  the  vertex  and  M  is  the 
mass  matrix. 

The  mass  matrix  arises  because  the  update  indicated  by  the  residual  R{W)  should  be 
made  to  the  average  value  in  the  control  volume.  It  thus  relates  the  average  value  of  a 
control  volume  associated  with  a  vertex  to  the  point  values  of  the  vertex  and  those  of  its 
immediate  neighbors.  This  definition  differs  from  the  way  the  mass  matrix  is  defined  in  finite 
element  formulations,  where  the  mass  matrix  arises  naturally  from  requiring  the  residual  of 
the  discretized  PDE  to  be  orthogonal  to  a  set  of  trial  functions,  with  the  solution  expanded  in 
a  set  of  basis  functions.  If  the  dissipative  terms  are  made  proportional  to  the  residuals,  this 
definition  of  the  mass  matrix  carries  through,  e.g.  the  Streamline  Upwind  Petrov  Galerkin 
method  [10].  For  finite  volume  schemes  employing  a  polynomial  reconstruction  procedure 
within  a  cell,  we  instead  derive  the  mass  matrix  entries  by  computing  the  average  of  this 
polynomial  over  the  control  volume.  The  mass  matrix  M  couples  the  system  of  ODE’s  in 
Eqn.  (3).  The  effect  is  that  even  with  an  explicit  scheme,  one  has  to  deal  with  the  solution 
of  a  coupled  linear  system.  A  technique  called  “mass-lumping”  [38],  replaces  the  matrix  M 
by  the  identity  matrix.  For  second-order  accurate  cell-centered  schemes,  which  employ  the 
triangles  as  the  control  volumes  and  store  the  values  at  the  centroids,  mass-lumping  does  not 
compromise  the  accuracy,  since  the  point  value  at  the  centroid  matches  the  average  value 
to  second  order.  However,  for  cell-vertex  schemes  on  nonuniform  grids,  the  centroid  of  the 
control  volume  is  not  represented  by  the  vertex  in  question.  For  time-accurate  computations 
on  such  grids,  mass  lumping  would  appear  to  introduce  locally  a  first  order  spatial  error. 
This  approximation  is  routinely  adopted  for  unsteady  flows  as  well,  and  does  not  appear  to 
adversely  affect  the  quality  of  the  solutions  obtained.  Davis  and  Bendiksen  [9]  observed  little 
discernible  differences  in  the  unsteady  solutions  when  using  the  full  and  the  lumped  mass 
matrices.  However,  since  they  used  an  explicit  scheme,  the  time  steps  were  quite  small  and 
furthermore,  the  grids  employed  appeared  to  be  fairly  uniform.  The  technique  employed  to 
solve  the  mass  matrix  (a  few  Jacobi  iterations)  in  [19,  9]  is  not  efficient,  especially  when  larger 
grids  are  used.  However,  Miller  [24]  and  Wathen  [45]  have  established  mesh-independent 
bounds  on  the  eigenvalues  of  the  diagonally  preconditioned  mass  matrix  arising  out  of  linear 
finite  elements  and  have  reported  that  a  conjugate  gradient  method  can  be  used  to  efficiently 
invert  the  diagonally  preconditioned  mass  matrix  in  just  a  few  iterations.  When  higher  order 
spatial  discretizations  are  employed,  the  mass  matrix  has  to  be  reckoned  with,  even  when 
using  cell-centered  discretizations.  We  note  that  the  mass  matrix  can  be  avoided  altogether 
if  only  cell  averages  are  employed  for  the  spatial  discretization. 

3  The  implicit  scheme 

We  first  outline  the  implicit  scheme  as  developed  by  Jameson  [11]  for  cell- centered,  structured 
grids,  where  mass  lumping  was  used.  Replacing  the  mass  matrix  in  Eqn.  (3)  by  the  identity 
matrix  and  making  a  3-point  backward-difference  approximation  for  the  time  derivative 
yields 

^  _ -I- 

2At  At  ^  2At 


3 


+R{W^+^)  =  0. 


(4) 


When  applied  to  a  linear  differential  equation  of  the  form, 

=  aW  (5) 

CLu 

this  discretization  is  A-stable  i.e.,  stable  for  all  values  of  a  At  in  the  left-half  of  the  complex 
plane  [8].  Eqn.  (4)  is  now  treated  as  a  steady  state  equation  by  introducing  a  pseudo-time 
variable  t*.  The  multigrid  scheme  then  solves  the  following  nonlinear  system  to  steady  state 
using  local  time  steps  At*: 

^  +  R-(U)  =  0,  (6) 

where  U  is  the  approximation  to  Here  the  unsteady  residual  R*[U)  is  defined  as 


R*{U)  =  ^VU  +  RiU)  - 


with  the  source  term 


n-lTJ/n-l 


=  — T/"W” - 


remaining  fixed  through  the  multigrid  procedure. 

A  multi-stage  Runge-Kutta  scheme  applied  to  solve  Eqn.  (6)  performs  the  role  of  a 
smoother.  A  low-storage  second  order  accurate  m-stage  Runge-Kutta  scheme  to  advance  U 
is  given  by 

go  =  u‘ 


y«+ig,  =  V^+^Qo-akAt*R*{Qk-i) 


=  Qm 

Starting  with  =  VE",  the  sequence  of  iterates  —  1,2,3....  converges  to 

However,  the  way  the  scheme  has  been  formulated  has  been  observed  by  Arnone  et  al.  [3] 
to  be  unstable  for  small  physical  time  steps.  At.  This  is  counter-intuitive  because  when 
using  a  small  At,  the  multigrid  procedure  should  converge  fast  and  ideally,  in  the  limit  of 
explicit  time  steps,  the  multigrid  procedure  should  converge  in  just  a  few  iterations.  Melson 
et  al.  [23]  showed  that  the  problem  is  due  an  instability  that  arises  when  a  small  At  is  used. 
They  modified  the  scheme  to  get  rid  of  this  instability.  The  source  of  the  problem  is  that 
the  unsteady  residual  R*(W)  includes  the  term  and,  is  therefore,  treated  explicitly 

in  the  Runge-Kutta  scheme.  Their  analysis  showed  that  if  this  term  were  treated  implicitly 
in  the  Runge-Kutta  scheme,  the  stability  region  would  grow  as  At  is  decreased.  It  is  easy  to 
treat  the  term  implicitly  since  it  is  only  a  diagonal  term.  Splitting  the  residual  R*(U)  as 

R"{u)  =  ^yu  +  R(U)  -  s,  (10) 
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the  Runge-Kutta  scheme  now  becomes, 


Qo  = 

I  +  (11) 

-akAr  [R{Qk-i)-S] 

[/'+'  =  Qm 

With  the  modified  scheme,  Melson  et  al.  [23]  have  shown  that  arbitrarily  large  or  small  At 
may  be  employed. 

As  in  [11,  23],  we  employ  a  full  approximation  storage  multigrid  scheme.  The  source 
term  is  computed  only  on  the  fine  grid  and  the  coarse  grid  problems  are  driven  by  the  fine 
grid  residuals.  For  the  generation  of  coarse  grids  we  follow  the  agglomeration  multigrid 
procedure.  In  this  method,  a  sequence  of  coarse  grids  is  generated  using  efficient  graph- 
based  algorithms.  This  method  has  certain  advantages  when  dealing  with  rigidly  moving  or 
deforming  meshes.  Since  the  edges  that  comprise  the  coarse  grid  volumes  are  subsets  of  the 
fine  grid  control  volume  edges,  when  the  grid  moves  rigidly  or  deforms,  the  projections  of 
the  control  volume  faces  onto  the  coordinate  directions  are  easily  computed  from  those  of 
the  fine  grid.  Also,  as  long  as  no  grid  points  are  added  or  removed,  and  the  triangulation 
remains  valid,  and  the  grid  connectivity  remains  unchanged,  the  interpolation  operators 
are  unchanged.  Multigrid  schemes  based  on  non-nested  triangulations  would  require  the 
recomputation  of  the  interpolation  operators  when  the  grids  deform.  Even  if  the  mesh 
connectivity  changes,  the  agglomeration  algorithm  is  efficient  enough  that  it  can  be  used  to 
regenerate  the  coarse  grids  without  incurring  substantial  overheads. 


4  Treatment  of  the  mass  matrix 


When  employing  a  vertex-centered  approximation,  making  a  3-point  backward-difference 
approximation  for  the  time  derivative  yields 

Q  2 

2At  At 


-f  =  0. 


The  multigrid  scheme  now  solves  the  following  system  to  steady  state  using  local  time  steps 

^  +  R*{U)  =  0,  (13) 

where  U  is  the  approximation  to  and  R*(W)  now  includes  the  mass  matrix  terms. 

Notice  that  the  first  term  -^VU  does  not  involve  the  mass  matrix,  thus  uncoupling  the 
system  of  equations.  The  explicit  Runge-Kutta  scheme  can  be  applied  exactly  as  before. 
The  inversion  of  the  mass  matrix  is  thus  accomplished  indirectly  during  the  multigrid  pro¬ 
cedure.  However,  the  modified  scheme  of  Melson  et  al.  [23]  poses  a  serious  problem.  Their 


modification  would  require  the  term  -^jVMU,  which  is  no  longer  a  diagonal  term,  to  be 
treated  implicitly.  We  have  devised  a  modification  that  solves  this  problem  which  is  detailed 
below.  The  implicit  Runge-Kutta  scheme  that  is  stable  for  all  At  is  given  by 

Qo= 

I  +  Qk  =  Qo  - 

a/iAt*  [R((5;fe-i)-5] 


=  Q 


(14) 


where  the  source  term  S  is  now  given  by 

S  =  - (15) 

At  2At 

If  we  simply  replace  the  mass  matrix  M  by  the  identity  on  the  left  hand  side  of  Eqn.  (14),  we 
have  observed  that  the  instability  at  small  time  steps  persists.  In  our  modification,  we  first 
add  and  subtract  ^^akAt*M^+^V^+^Qk-i  on  the  right  hand  side  of  Eqn.  (14)  to  obtain 


I  +  -^akAt*M^+A  V^+^Qk  =  V^+'^Qo  - 
akAt*R*{Qk-i)  +  ^akAt*M^+^V-^^Qk-i,  (16) 


where  use  has  been  made  of  the  equation 

r*(U)  =  ^VMU  +  R{U)-S.  (17) 

Note  that  the  same  term  -^^ockAt* MV Q  appears  on  the  left  and  the  right  hand  sides  of 
Eqn.  (16),  except  that  they  are  evaluated  at  the  h  and  k  stages,  and  that  R*{U)  is  being 
driven  to  zero.  The  mass  matrix  M  can  now  be  replaced  by  where  I  is  the  identity 
matrix  and  /?  is  a  constant  yielding  the  following  equation: 

-a,Af  R"{Q„-,)  +  (18) 

The  method  can  always  be  stabilized  by  increasing  /?  and  is  akin  to  using  a  damped  Jacobi 
method.  The  Runge-Kutta  scheme  no  longer  requires  a  matrix  inversion.  For  small  time 
steps  of  the  order  permitted  by  the  explicit  scheme,  we  find  that  the  choice  of  /?  =  2  stabilizes 
the  scheme. 
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5  Mesh  point  movement  strategies 

In  order  to  be  able  to  perform  unsteady  flow  simulations  over  moving  geometries,  a  body- 
conforming  mesh  has  to  be  regenerated  either  globally  at  each  time  step,  or  the  existing  grid 
can  be  allowed  to  deform.  The  former  option  is  expensive,  especially  in  three  dimensions. 
However,  Baum  et  al.  [6]  have  proposed  some  simplifying  strategies.  Whenever  the  body 
movement  becomes  too  severe,  they  regenerate  a  coarse  mesh  either  locally  or  globally.  This 
is  followed  by  the  use  of  adaptive  h-refinement  in  order  to  create  a  fine  mesh.  In  the  present 
work,  we  only  investigate  mesh  point  movement  strategies.  Tension  spring  analogy  [5,  28,  35] 
or  other  physical  analogies,  such  as  incompressible  flow  [12],  are  typically  used  to  move  the 
mesh  points.  In  the  former  case,  the  distribution  of  the  spring  stiffnesses  is  crucial.  Since  the 
techniques  try  to  maintain  the  connectivity  of  the  grid  at  all  costs,  the  grid  lines  may  cross 
resulting  in  invalid  triangulations.  Nevertheless,  for  many  simple  configurations,  such  as 
isolated  airfoils,  no  crossover  occurs  and  the  spring  analogy  has  been  demonstrated  to  work 
well  e.g.,  [5].  For  such  cases  however,  the  use  of  exponentially  varying  scaling  factors  [9]  or 
rigid  body  motion  are  simpler.  The  real  challenge  for  any  mesh  point  movement  strategy  is 
when  relative  motion  is  present  between  bodies  in  close  proximity  and  when  fine  grids  are 
employed.  Other  possibilities  for  grid  movement  are  methods  in  use  in  the  moving  finite 
element  method  [25],  where  evolution  equations  are  derived  for  grid  point  motion  from  the 
governing  equations. 

In  our  approach,  we  use  the  tension  spring  analogy  to  allow  the  grid  points  to  react  to 
the  motion  of  the  geometries.  The  following  linear  system  of  equations  is  solved  by  a  Jacobi 
method: 

^  Kij{xi  —  Xj)  =  Sij,  (19) 

3 

where  the  source  term  Sij  is  computed  so  as  to  maintain  the  initial  grid  in  the  absence  of  any 
displacements.  The  spring  stiffness  Kij  is  taken  to  be  where  lij  is  the  length  of  the  edge 
joining  nodes  i  and  j.  Over  a  number  of  applications,  we  have  found  the  value  of  p  =  2  to 
work  well.  When  the  boundary  motion  is  prescribed,  this  method  does  not  guarantee  that 
the  grid  lines  will  not  cross.  The  next  improvement  is  to  make  the  spring  system  nonlinear 
i.e.,  the  boundary  motion  is  decomposed  into  smaller  steps  and  the  procedure  is  repeated  at 
every  step.  When  relative  motion  is  present,  this  results  in  excessive  skewing  of  grid  lines  and 
eventually  the  lines  do  cross.  Thus  some  kind  of  reconnection  procedure  is  unavoidable.  We 
choose  to  swap  the  edges  by  the  Delaunay  criterion.  However,  even  though  the  grid  remains 
valid,  the  resulting  point  distribution  is  typically  unsatisfactory.  Therefore  a  presmoothing 
procedure  is  used  to  smooth  the  distribution  of  points.  This  involves  using  a  Jacobi  method 
on  the  following  system: 


(/  +  =  xf  -F  (20) 

3 

where  Ni  is  the  degree  of  node  i.  Typically  we  find  that  of  the  order  of  100-200  steps 
are  required  to  solve  Eqn.  (19),  while  only  about  3-4  iterations  of  Eqn.  (20)  are  sufficient. 
The  number  of  Jacobi  iterations  to  solve  Eqn.  (19)  may  appear  to  be  large,  but  this  is 
mainly  due  to  the  large  displacements  arising  from  using  large  time  steps  permitted  by  the 
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implicit  scheme.  The  grid  motion  procedure  described  above  also  has  applications  in  design 
optimization,  where  the  surface  geometry  changes  during  the  design  cycle. 

For  the  implicit  flow  solver,  the  coarse  grids  are  rederived  by  using  the  agglomeration 
algorithm  whenever  the  grid  connectivity  changes.  An  incremental  agglomeration  algorithm 
is  possible  that  operates  only  on  the  affected  region,  but  this  is  not  currently  done.  Typically, 
for  the  time  steps  used  with  the  implicit  scheme,  the  grid  reconnection  and  reagglomeration 
consume  about  a  third  of  the  time  required  for  the  flow  solver,  mainly  on  account  of  the 
swapping  procedure  not  being  vectorized. 

The  grid  movement  terms  need  to  be  discretized  carefully  so  that  freestream  is  preserved. 
In  other  words,  simply  moving  the  grid  through  the  domain  should  not  change  the  freestream 
solution.  The  Geometric  Conservation  Law  (GCL)  [39,  47,  27]  formalizes  this  concept.  It 
can  be  derived  from  the  continuity  equation  in  Eqn.  (1)  by  first  assuming  the  control  volumes 
to  be  the  simplices  themselves.  Assuming  a  uniform  velocity  field  and  a  constant  density 
field,  we  obtain 

^+<f  [V  —  s].n  da  =  0,  (21) 

at  Js{t) 

where  V  is  the  velocity  field  and  s  is  the  velocity  of  the  boundary  S{t).  Since  V  is  constant 
and  the  control  volume  is  assumed  to  be  closed  at  all  times  so  that  n  da  =  0,  the 
equation  becomes 

—  i  s.n  da  =  0.  (22) 

dt  Js(t) 

The  discrete  form  of  this  equation  should  hold  at  all  time  steps  and  for  all  the  simplices  and 
is  called  the  GCL.  Using  a  forward  Euler  approximation  for  the  time  derivative,  we  obtain 


rt+At  /• 

_  V”  =  /  (f)  s.n  da  dt 

Jt  Jsi(t) 

_  /‘t+At  r 


/‘t-f-m  r 

=  y  /  s.n  da  dt., 

:  Jt  JUj  At) 


where  Sj{t)  =  Ylj  ^/,i(^)  is  the  surface  enclosing  the  volume  V/(t)  of  simplex  I.  The  term 
inside  the  summation  represents  the  volume  swept  out  by  the  boundary  fl/j  as  the  grid 
points  forming  that  segment  move.  If  the  grid  points  are  allowed  to  move  arbitrarily,  the 
GCL  enables  the  edge  velocity  s  and  the  grid  normal  n  to  be  determined.  Since  simplices 
are  convex,  the  volumes  V”,  are  uniquely  determined  by  the  positions  of  the  points  at 
time  levels  n  and  n  +  1.  In  two  dimensions,  it  can  be  shown  [47]  that  the  change  in  volume 
can  be  expressed  as 

^n+l  _  yn  =  At  ^  (24) 


where  the  summation  is  over  the  edges  forming  the  triangle  /.  Here  the  edge  velocity 
is  given  by  the  average  of  the  velocities  of  the  vertices  connected  by  the  edge  and  N the 
normal  scaled  by  the  edge  length,  is  given  by  the  average  of  the  values  at  the  old  and  new 
time  level: 

=  0.5(Arj  +  iV”+'). 


(25) 


The  velocity  is  assumed  to  be  constant  within  a  time  step.  Therefore,  the  velocity  of  vertex 
i  is  given  by 

Y-n+l 

\  (26) 

where  X  is  the  position  vector  of  the  vertex. 


Figure  1:  Control  volume  for  vertex  i  and  edge  normals. 


For  a  vertex-centered  scheme,  such  ais  the  one  used  in  the  present  work,  the  control 
volumes  are  formed  by  the  median  dual.  We  closely  follow  the  excellent  development  of 
Nkonga  and  Guillard  [27],  who  have  presented  the  steps  involved  in  properly  discretizing  the 
flow  equations  in  three  dimensions  so  as  to  obey  the  GCL  for  a  two-level  scheme  in  time.  The 
control  volume  edges  are  delimited  by  segments  of  the  medians  of  the  triangles  (Figure  1). 
Inspecting  one  edge  of  the  triangle,  we  notice  that  the  control  volume  edge  is  comprised  of 
two,  generally  non-collinear  median  segments.  The  results  from  Eqns.(24,  25)  can  be  used  to 
express  the  change  in  the  volume  of  any  polygon  cis  a  summation  over  the  edges  comprising 
the  polygon  of  the  volumes  swept  out.  We  obtain  the  following  expression  for  the  change  in 
Vi,  the  volume  associated  with  vertex  i: 


yn+l  _  yn  ^  At  ^ 


Tl-f-l/2  71+1/2 


where  the  summation  is  over  the  edges  meeting  at  vertex  i.  Also,  and  A/'JJ 


normals  scaled  by  the  lengths  to  the  left  and  right  of  the  edge,  are  obtained  as  averages  of 
the  values  at  n  and  n-|- 1  time  levels.  The  edge  velocities,  and  are  computed  as 

the  averages  of  the  velocity  of  the  midpoint  of  the  edge  and  that  of  the  left  or  right  centroid, 
respectively. 

In  a  vertex-centered  finite- volume  setting,  the  discretization  of  the  governing  equations 
(Eqn.  (1))  for  a  two-level  explicit  or  implicit  scheme  that  obeys  the  GCL  is  performed  as 
follows: 

{Mvwy^^  -  {Mvwy  ^  ^  _  ^n+1/2^  _  5«+i/2y 


T  71+1/2 


)  +  K 


n+1/2 
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where  t)  =  n  for  forward  and  rj  =  n-\-l  for  backward  Euler  methods,  and  W  is  the  average  of 
the  values  at  the  vertices  connected  by  the  edge  j ,  Also,  gj  represents  the  discretization  of 
the  last  term  in  Eqn.  (1)  and  makes  use  of  the  scaled  normal  for  the  edge,  ^  ,  defined 


^n+l/2  _  ^n+1/2 


When  a  three  point  backward  difference  approximation  is  used  for  the  time  derivative, 
the  spatial  discretization  is  performed  as  follows: 

Z{MVWY^^  -4{MVW)^  +  [MVWf-^ 

2At 

+  -  ‘’it''")  +  -  ‘Tk'")) 


•n-1/2  /-rrn+l 


=  0. 


It  is  easy  to  see  that  this  formulation  satisfies  the  GCL,  and  that  it  reduces  to  the  standard 
three-point  formula  (Eqn.  (12))  in  the  absence  of  any  grid  motion.  We  mention  here  that 
for  fixed  grid  applications  it  is  not  necessary  to  compute  N and  N j^R  separately.  Instead 
only  the  sum  of  the  two  is  required,  which  can  be  computed  easily  from  the  the  centroidal 
coordinates.  On  the  other  hand,  for  moving  grid  problems,  the  N and  N j^R  need  to  be 
computed  separately  for  use  in  Eqns.  (28,30). 

The  spatial  discretization  given  above  only  deals  with  the  convective  and  diffusive  fluxes. 
The  dissipative  fluxes  are  computed  in  the  usual  fashion  using  the  scaled  normal  given  by 
Eqn.  (29)  for  the  two-level  explicit  or  implicit  scheme.  For  the  three-point  scheme,  the  scaled 
normal  is  defined  as  follows: 

TV,-  =  3/2  -I-  -  1/2  .  (31) 


TV,  =  3/2  -  1/2  .  (31) 

The  discretization  of  the  dissipative  fluxes  has  no  implications  for  the  GCL  as  these  fluxes 
vanish  for  a  uniform  field. 

We  have  verified  that  the  formulation  given  above  preserves  the  freestream  conditions 
to  machine  precision.  With  other  treatments,  such  as  simply  modifying  the  fluxes  at  the 
vertices  by  using  nodal  velocities,  freestream  conditions  are  preserved  to  much  less  precision, 
which  could  be  detrimental  in  some  applications. 

Finally,  we  address  the  issue  of  boundary  conditions  for  inviscid  flows.  The  control  volume 
for  a  boundary  vertex  is  illustrated  in  Figure  2.  The  summation  over  the  edges  in  Eqns. 
(28,30)  is  augmented  with  the  contributions  from  the  the  boundary  edges  that  delimit  the 
control  volume.  For  the  two- level  scheme,  the  following  terms  are  added  to  Eqn.  (28): 


w2iNn' 


.n+l/2 


)i + '')] 


+  (32) 

where  Wi  and  Wr  are  the  values  of  W  on  the  left  and  right  side  of  vertex  i,  Nh,R  and  Nb,L 
are  the  outward  normals  to  the  boundary  edges  of  the  control  volumes,  and  and  s^^r 
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b,L 


Figure  2:  Control  volume  for  vertex  boundary  vertex  i  and  boundary  edge  normals. 


are  the  grid  velocities  associated  with  these  edges.  For  an  inviscid  all,  since  the  [V  —  s]  =  0, 
only  the  last  two  terms  are  retained.  Likewise  for  the  three-point  scheme  (Eqn.  (30)),  the 
additional  terms  are 


+  '')i 

-  (33) 

For  an  inviscid  wall,  only  the  last  four  terms  are  retained,  implying  the  following  boundary 
condition: 


-jv 


_  3  ^n+l/2  n+1/2  1  i^n-ll2  n-1/2 


6  Mesh  restructuring  and  interpolation 

The  procedure  for  mesh  restructuring  was  discussed  earlier.  This  involved  a  mechanism  for 
moving  grid  points  in  response  to  the  motion  of  the  boundary  points,  evolving  the  solution 
with  the  grid  movement  terms  present,  and  swapping  the  edges  of  the  triangulation  to 
improve  grid  quality.  When  using  the  three  point  difference  formula  (Eqn. (12)),  the  swapping 
of  edges  at  a  particular  time  step  has  to  be  done  in  such  a  manner  that  the  triangulation  is  also 
valid  (i.e.,  no  crossing  of  edges)  at  the  previous  time  step.  This  additional  criterion  can  easily 
be  incorporated  in  the  swapping  procedure.  Using  the  trapezoidal  rule  would  circumvent 
this  step  and  still  ensure  second  order  accuracy  in  time.  However  it  is  unattractive,  because 
of  the  instabilities  that  occur  as  large  time  step  sizes  are  used  [8]. 
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When  the  swapping  of  edges  occurs,  the  solution,  which  is  stored  at  the  vertices  of  the 
triangulation,  needs  to  be  modified.  For  capturing  weak  solutions  of  hyperbolic  conserva¬ 
tion  laws,  conservation  in  time  is  a  requirement.  If  conservation  were  not  an  issue,  the 
solution  need  not  be  modified  since  it  merely  changes  from  one  piecewise-linear  representa¬ 
tion  to  another  and  the  two  are  acceptable  second-order  accurate  solutions.  Note  that  the 
swapping  takes  place  at  a  fixed  time  step  after  the  solution  has  already  been  evolved  up 
to  that  time  level.  It  is  possible  to  carry  out  another  iterative  loop  so  that  the  unsteady 
residual  is  driven  to  zero  in  the  new  configuration.  This  would  double  the  work  done  at 
each  time  step.  Instead,  we  propose  a  noniterative  conservative  interpolation  at  each  time 
step.  With  the  lumped  mass  matrix,  Eqn.  (4),  the  requirement  for  conservation  is  that 
Uvertices  ^0  conserved  before  and  after  the  swapping.  When  an  edge  is  swapped,  the 
control  volumes  change  for  all  the  four  points  forming  the  quadrilateral.  Figure  15  shows 
a  four  point  quadrilateral  subset  of  a  triangulation.  The  initial  triangulation  given  by  the 
triangles  125  and  126  and  the  triangulation  after  swapping  given  by  the  triangles  156  and 
256.  The  portions  of  the  control  volumes  associated  with  the  four  vertices  that  fall  inside  the 
quadrilateral  region  are  illustrated  as  well.  A  conservative  interpolation  can  be  performed 
by  treating  the  solution  inside  the  control  volumes  as  piecewise  constants  and  computing 
geometrically  the  fractions  of  the  old  control  volumes  comprising  the  new  ones.  This  would 
require  complex  intersection  computations,  especially  for  the  vertex-centered  scheme,  where 
the  control  volume  edges  are  segmented  edges.  This  algorithm  may  be  viewed  as  a  particular 
application  of  the  algorithm  due  to  Ramshaw  [32].  An  algebraic  approach  that  only  involves 
areas  Ai,  A2,  A5,  Aq,  Ay^  Ay,  A5/,  A^t  is  attractive,  but  the  system  becomes  underdetermined 
if  conservation  is  assumed.  The  drawbacks  of  these  algorithms  are  algorithmic  complexity 
and  their  diffusive  character.  The  latter  is  particularly  nettlesome  because  it  degrades  the 
accuracy  to  first  order.  For  example,  swapping  back  to  the  initial  configuration  and  repeating 
the  procedure  in  Figure  15,  will  result  in  equal  values  at  the  four  vertices. 


5  5 


6  6 


OLD  NEW 

Figure  3:  Two  possible  triangulations  of  four  points  and  control  volumes  of  vertices. 


12 


In  the  following,  a  conservative,  linearity-preserving  interpolation  procedure  is  presented. 
The  derivation  is  given  assuming  that  the  lumped  mass  approximation  is  adopted.  We 
observe  that  the  following  identity  holds: 

^  W= 

Vertices  Triangles 

where  V  is  the  area  of  the  control  volume  and  V  is  the  area  of  the  triangle,  and  W,  the 
average  value  in  a  triangle,  is  given  by  the  average  of  the  vertex  values.  Referring  to  Figure 
15,  when  an  edge  is  swapped,  the  contributions  from  the  triangles,  125  and  126,  change  to 
the  contributions  from  the  two  swapped  triangles,  156  and  256,  on  the  right  hand  side  of 
Eqn.  (35).  We  further  note  that  each  term  in  the  right  hand-side  may  be  viewed  as  volume 
under  a  “roof”  inthe  x  —  y  —  W  space.  If  the  data  were  linear,  the  volumes  under  the  two 
triangulations  would  be  identical.  In  that  case,  no  changes  need  to  be  made  to  the  solution 
at  the  vertices.  In  the  general  case,  we  compute  the  contributions  to  the  right-hand  side  of 
Eqn.  (35)  from  the  two  triangulations  as: 

Tm  =  ^[{Wi  +  W2  +  Ws)Vr25  +  iWi+W2  +  We)Vi2e]. 

Tnev.  =  \[iWr  +  Ws  +  We)Vise  +  iW2  +  Ws  +  We)V2B6].  (36) 

Changes  are  now  made  to  the  vertex  values  in  a  conservative  manner  by  distributing  the 
difference.  Told  —  Tnew,  to  the  nodes.  In  principle,  the  changes  can  be  made  to  all  the  four 
vertex  values  and  would  still  result  in  a  linearity-preserving,  conservative  scheme,  but  in 
inspecting  Figure  15  we  notice  that  with  swapping,  the  control  volumes  associated  with 
vertices  1  and  2  (connected  by  the  old  edge)  can  only  decrease,  whereas  those  associated 
with  vertices  5  and  6  (connected  by  the  new  edge)  can  only  increase.  Therefore  changes  need 
to  be  made  only  to  vertices  5  and  6.  We  apportion  the  difference,  Toid  —  Tnewi  equally: 

^  -f  i(T<,W  ”  T^e^), 

+  ^{Told  -  Tne^)-  (37) 

We  now  show  that  the  interpolation  formulas  given  by  Eqn.  (37)  satisfy  the  conservation 
property.  First  note  that  the  overlapping  control  volume  given  by  the  union  of  the 
triangles  meeting  at  vertex  i,  is  related  to  the  nonoverlapping  control  volume  Vi  by  the 
formula, 

vi=^3Vi  =  J2V^  (38) 

j 

where  the  sum  is  over  all  the  triangles  meeting  at  vertex  i.  After  swapping,  we  have 

=  ^5  -  Vm  +  Vise  +  V256  =  ^5  +  Vi26.  (39) 

We  can  thus  derive  the  following  relations: 

yneu,  ^  V,  +  Vus/^. 

y--  =  V6  +  V125/3, 

ynevj  ^  ^  -  V256/3, 

yne^  =  ^2  "  ViseA 


With  the  help  of  these  relations,  it  is  easy  to  show  that  the  conservation  property  holds: 

^^newynew  ^  yy^ynen,  ^  yy^y^  y  yy^y^  y  yy^y^^  (40) 

where  use  hcis  been  made  of  the  fact  that  the  values  at  vertices  1  and  2  do  not  change. 

In  the  formulas  given  by  Eqn.  (37),  the  difference  Toid—Tnew  has  been  apportioned  equally 
to  the  two  vertices.  It  is  easy  to  see  that  this  could  introduce  new  extrema  in  the  solution. 
The  degree  of  freedom  available  in  how  this  apportionment  is  made  to  the  two  vertices  can 
be  used  to  effectively  prevent  new  extrema.  Let 


^^min  — 


=  Max(Wi,W2,W5,W6) 
=  Min(Wi,W2,W5,W6). 


For  e  =  5  or  6  compute. 


\iToli-Tner.>^ 

~  1  ynew  [{  T  U  -  T  <0 

I  Toid-Tncw'  -Lnew<~yf- 

Let  r  =  Min[0.5,r5,r6].  The  interpolation  formulas  can  be  written  in  a  compact  form  as 
follows: 


lyneu) 


=  W5  + 


T'old — T'new 

ynetu 


{r-rs}{r-re}  n  c  ,  (r-U.5)(r-r6 J  ,  / 

(0.5— r5)(0.5-r6)  *  ‘  (rs— 0.5)(r5— re)  '  (re— 0.5)(r6—r5) ' 

(r-r5)(r-r6)  a  c  i  (r-0.5)(r-r6)  /-■  _  ..'j  1 
(0.5— r5)(0.5— re)  '  '  (rs— 0.5)(r5— re)  '  '  (’'6— 0.5)( 


Ti/neii)  _  Tir  I  Tm-Tu^w  \  (r-rs)(r-re)  a  c  i  (r-0.5)(r-re)  /-■  \  ,  (r-0.5)(r-r5j_  J 

W/g  -  H/g  +  Tneu,  [(O.S-rsKO.S-re)^-^  +  (r5-0.5)(r5-re)  +  (re -0.5) (re -re )n  ’ 

The  formulas  are  conservative,  linearity-preserving,  and  also  do  not  introduce  new  extrema. 
In  contrast  to  the  conservative  formulas  that  assume  piecewise  constant  values  and  result  in 
equal  values  at  the  vertices  upon  repeated  application,  the  new  formulas  converge  locally  to 
a  linear  (possibly  least  squares)  profile  for  the  four  data  points  upon  repeated  application. 
The  swapping  and  the  interpolation  formulas  generalize  to  three-dimensional  tetrahedral 
tessellations  and  are  discussed  in  the  Appendix. 


7  Results 

First,  results  from  a  one- dimensional  example  are  presented  illustrating  the  role  of  the  mass 
matrix.  On  a  uniform  mesh,  since  the  vertex  and  the  centroid  of  its  control  volume  coincide, 
the  mass  matrix  can  be  lumped  without  suffering  any  adverse  consequences.  The  situation  is 
different  if  a  mesh  with  variable  mesh  widths  is  considered.  The  one-dimensional  advection 
equation 

¥  +  =  («> 


is  solved  on  a  random  grid.  The  scheme  stores  the  pointwise  values  u*  at  locations  Xi.  The 
initial  condition  is  a  Gaussian  and  the  profile  is  advected  by  marching  to  a  fixed  time.  A  grid 
refinement  study  is  carried  out  using  a  constant  CFL  number  of  0.3.  The  spatial  derivative 
is  approximated  in  a  MUSCL  scheme  [40]  as 
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(43) 


where  the  value  interpolated  to  the  left  side  of  the  face  i  +  1/2,  is  given  by 


^i+l /2 


Ui  +  (a^i+1/2  ~ 


'^i—1 


The  mass  matrix,  which  is  tridiagonal,  is  inverted  using  the  Thomas  algorithm.  We  have 
experimented  with  two  definitions  of  the  mass  matrix.  The  first  one  derives  the  mass  matrix 
by  assuming  a  piecewise  linear  distribution  of  data  between  the  grid  points  and  computes 
the  average  over  the  control  volume  a:^+i/2]- 


1  3  1  Xi^i  Xx 

- Ui^l  +  -Ui  +  - - 

Xi—i  4  4  x^^i 


A  second  definition  of  the  mass  matrix  is  derived  by  computing  the  average  of  the  recon¬ 
struction  polynomial  within  a  control  volume.  This  polynomial  is  given  by 

u{x)  =  Ui  +  (a:  -  (46) 


The  average  over  the  control  volume  is  given  by 

Xx—'i  4"  ,  Xi^i  { A7\ 

TT  '  r  1  “h  T  .  N  y^i—i  V  / 

8(xi4.i  +  Xi-i)  8(x,-+i  +  Xi^i) 

Figure  1  compares  the  errors  in  L2  norm  with  the  mass  matrices  given  by  Eqn.  (45)  and 
Eqn.  (47),  and  with  the  lumped  mass  matrix.  All  the  schemes  exhibit  second  order  accuracy 
but  the  errors  are  larger  with  the  mass  matrix  given  by  Eqn.  (45).  The  results  obtained 
with  the  lumped  mass  matrix  are  almost  identical  to  those  obtained  with  Eqn.  (47).  Taylor 
series  expansion  would  imply  a  first  order  error  with  the  lumped  mass  matrix  on  a  random 
grid,  whereas  Figure  1  clearly  indicates  second  order  accuracy.  The  results  therefore  reveal 
the  inadequacy  of  local  analysis.  The  results  obtained  with  the  usual  finite  element  mass 
matrix,  with  hi  =  Xi  —  Xi_i, 

2/i^  4  ^  2hi^i  ^ 

— - — — - rUi-i  +  -Ui  +  —7  --7  j 

6[hi  +  hi^i)  6  6{hi  +  hiA-i) 

are  also  shown  in  Figure  1  and  again  display  larger  errors  compared  to  the  lumped  mass 
approximation.  The  reason  for  this  is  that  the  finite  element  mass  matrix  is  consistent 
with  a  Galerkin  method  which  corresponds  to  a  central  difference  discretization,  whereas 
the  spatial  differencing  employed  here  is  upwind-biased.  After  experimenting  with  a  one- 
parameter  family  of  mass  matrices,  we  have  found  that  the  lumped  mass  matrix  gives  the 
lowest  errors  with  this  particular  spatial  discretization. 

It  is  well  known  in  finite  element  literature  [38]  that  in  some  cases  the  lumping  of  the 
mass  matrix  does  not  compromise  the  solution  accuracy,  but  that  the  mass  matrix  may  play 
a  crucial  role  when  higher-order  discretizations  are  considered.  To  examine  this,  we  employ 
a  quadratic  reconstruction  procedure  utilizing  point  values.  With  hi  =  Xi  —  Xi^i^  we  obtain. 


7/^ 
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^hi{hi  +  hij^i) 
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+  TTT - 
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Figure  4:  L2  norms  of  the  errors  with  various  schemes  on  a  random  grid 


gure  5:  L2  norms  of  the  errors  with  quadratic  reconstruction  on  a  random  grid. 


The  finite  volume  mass  matrix  is  derived  by  determining  the  average  of  the  quadratic  dis¬ 
tribution  and  is  given  by: 

aiUi-\  -)-  a2Ui  oaUj-i-i, 


2hihi^i  -|-  Kf 
12hi[hi 

12hi^i{hi  +  hi+i) 

ct2  —  1  —  —  ®3‘ 


The  standard  Runge-Kutta  scheme  which  is  fourth  order  accurate  in  time  is  used  for  the 
higher  order  computations.  Figure  2  shows  the  error  plots  using  the  lumped  mass  matrix 
and  the  full  mass  matrix  on  the  random  grid.  It  shows  that  with  the  lumped  mass  matrix, 
the  accuracy  eventually  degrades  to  second  order  as  the  grid  is  refined,  whereas  using  the 
full  matrix  yields  the  third  order  accuracy  of  the  spatial  discretization.  We  have  observed 
that  using  any  other  definition  for  the  mass  matrix  degrades  the  accuracy  to  second  order. 

The  implications  for  the  scheme  in  multiple  dimensions  are  clear.  As  long  as  only  a 
second  order  accurate  scheme  is  used  and  we  operate  with  either  cell-vertex  or  cell-centered 
data,  the  mass  matrix  may  be  lumped  without  any  loss  of  order  of  accuracy.  The  mass 
matrix  can  also  be  ignored  for  second  (and  higher)  order  accurate  schemes  if  a  strict  cell- 
average  interpretation  is  employed.  If  point  values  are  used  to  construct  third  and  higher 
order  accurate  schemes,  the  accuracy  will  degrade  if  the  mass  matrix  is  lumped.  For  higher 
order  accurate  schemes  based  on  point  values,  the  indirect  mass  matrix  inversion  technique 
discussed  earlier  will  help  preserve  the  order  of  accuracy  of  the  scheme. 

We  next  present  results  from  two-dimensional  inviscid  calculation  over  a  pitching  airfoil. 
The  transonic  flow  is  over  a  sinusoidally  oscillating  NACA0012  airfoil  where  the  angle  of 
attack  a{t)  varies  according  to  the  formula 

a{t)  =  Qtm  +  Ciosin{u)t)  (49) 

For  the  test  case  chosen,  am  =  0.016°,  ao  =  2.51°,  k  =  ^  =  0.0814  and  the  freestream 
Mach  number.  Moo  =  0.755.  Computing  this  flow  using  an  explicit  scheme  is  time-consuming 
because  of  the  low  frequency.  The  flows  are  computed  using  two  meshes,  referred  to  as  GRIDl 
and  GRID2,  each  having  6336  vertices.  These  are  shown  in  Figures  3  and  4,  respectively. 
GRIDl  is  generated  by  drawing  diagonals  in  a  structured  C-mesh  and  is  fairly  uniform. 
GRID2  is  generated  by  random  perturbations  on  GRIDl.  In  GRID2,  the  centroids  of  the 
control  volumes  formed  by  the  median  dual  are  not  represented  well  by  the  vertices.  Figure  5 
shows  the  lift  histories  during  the  third  cycle  of  oscillation.  Four  curves  are  shown,  namely, 
the  histories  with  the  lumped  and  full  mass  matrices  for  GRIDl  and  GRID2.  Since  a 
projection-evolution  approach  is  not  followed,  the  mass  matrix  is  derived  by  using  a  definition 
similar  to  Eqn.  (45).  As  expected,  the  mass  matrix  has  little  impact  on  the  integrated 
quantities  even  in  the  random  mesh.  The  differences  in  the  solutions  between  the  two  grids 
are  likewise  insignificant.  The  CPU  time  increases  by  about  15%  when  the  full  mass  matrix 
is  included.  These  examples  have  been  run  with  a  maximum  physical  CFL  number  of  500, 
corresponding  to  using  54  time  steps  per  sinusoidal  oscillation  of  the  airfoil.  The  number  of 
iterations  for  the  inner  multigrid  procedure  is  fixed  at  30.  Figure  6  shows  the  convergence 
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of  4-grid  agglomeration  multigrid  procedure  during  a  particular  time  step  with  the  lumped 
and  the  full  matrices  where  the  L2  norm  of  the  unsteady  residual  R*  is  plotted  as  a  function 
of  the  multigrid  cycles  using  4  grid  levels.  The  convergence  improves  slightly  when  the  mass 
matrix  is  included.  The  reason  for  this  is  that  the  mass  matrix  can  be  recast  as  an  implicit 
smoothing  operator.  The  remaining  examples  have  been  computed  with  the  lumped  mass 
matrix.  Finally,  Figures  7  and  8  shows  the  effect  of  the  physical  time  step  size.  Three 
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Figure  8:  Lift  histories  during  the  third  cycle  of  motion. 

lift  and  moment  histories  are  shown,  employing  40,  20  and  10  steps  per  pitch  cycle  using 
the  lumped  mass  matrix  and  deforming  grids.  The  grid  velocities  are  computed  by  using 
Eqn.  (27).  The  nodes  are  repositioned  by  the  procedure  described  earlier.  No  swapping 
of  edges  occurred  in  this  computation.  The  number  of  multigrid  cycles  used  at  each  time 
step  are  respectively,  15,  20  and  20.  With  10  steps  per  pitch  cycle,  some  discrepancy  may 
be  observed;  increasing  the  number  of  multigrid  cycles  (thus  solving  the  nonlinear  problem 
better  at  each  time  step)  did  improve  the  solution.  Also  shown  is  a  comparison  with  the 
experimental  data  of  Landon  [16]  as  well  as  with  a  structured  grid  computation  employing 
a  TVD  scheme  on  the  same  grid  [42].  The  differences  in  the  moment  histories  between  the 
structured  and  unstructured  grid  computations  may  be  attributed  to  the  differences  in  the 
spatial  discretization.  To  assess  the  effects  of  swapping,  we  modify  the  Delaunay  criterion  to 
force  the  swapping  of  edges.  Given  a  quadrilateral  with  two  possible  triangulations,  rather 
than  swap  to  the  new  configuration  if  it  has  a  larger  minimum  angle,  we  swap  if  some 
multiple  (1.1  in  this  example)  of  the  minimum  angle  in  the  new  configuration  is  larger  than 
the  minimum  angle  in  the  old  configuration;  only  one  pass  of  this  algorithm  is  performed. 
On  a  coarse  triangular  grid  generated  from  a  structured  128  x  32  grid,  the  moment  histories 
during  the  first  three  pitch  cycles  without  the  swapping  are  shown  in  Figure  9.  Figure 
10  shows  the  histories  during  the  first  five  pitch  cycles  with  the  swapping  of  edges.  It  is 
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seen  that  the  transient  response  is  indeed  quite  different  especially  during  the  second  cycle, 
although  the  response  does  become  periodic  during  the  fifth  cycle.  Figure  11  shows  the 
moment  histories  on  a  triangular  mesh  generated  from  a  256  x  64  grid  with  the  swapping  of 
edges.  The  transient  response  is  not  nearly  as  erratic  as  on  the  128  x  32  grid  with  swapping 
and  more  in  line  with  the  response  on  the  coarse  grid  without  swapping.  Note  that  the 
swapping  and  subsequent  conservative  interpolation  introduces  errors  proportional  to  Aa;^, 
which  decrease  as  the  grid  is  refined. 

We  next  present  a  laminar  unsteady  calculation  of  impulsive  start-up  flow  over  a  cylinder 
at  a  Reynolds  number  of  1200.  Rumsey  [34]  has  computed  this  flow  using  a  structured  grid 
code  and  has  made  detailed  comparisons  with  experimental  data  [26].  We  employ  the  same 
grid  (192  X  64)  as  in  [34],  but  divide  the  quadrilaterals  into  triangles.  The  flow  separates 
and  eventually  vortices  are  shed.  The  computed  Strouhal  number  is  0.225  as  compared  to 
the  experimental  value  of  0.215  and  the  computed  value  of  0.222  using  the  structured  grid. 
The  time  is  nondimensionalized  as  t  =  t/(d/f/oo)  where  d  is  the  diameter  of  the  cylinder 
and  Uoo  is  the  freestream  velocity.  The  sequence  of  non-dimensional  time  steps  At  is  chosen 
as  in  [34]:  At  =  0.01  for  t  up  to  6.0,  At  =  0.02  between  t  =  6.0  and  9.0  and  At  =  0.05 
above  t  =  9.0.  We  use  the  agglomeration  multigrid  strategy  with  six  grids  where  the  edge- 
coefficients  needed  for  computing  the  viscous  and  inviscid  terms  on  the  coarse  grids  are 
precomputed  as  in  [22].  Seven  multigrid  iterations  were  used  at  each  time  step  yielding 
about  of  2-3  orders  of  reduction  in  the  unsteady  residual.  The  centerline  velocity  at  a  time 
t  =  2.9  is  plotted  in  Figure  8,  and  the  maximum  reverse  flow  velocity  is  plotted  in  Figure  9 
as  a  function  of  i  during  the  initial  bubble  growth  phase.  Good  agreement  may  be  observed 
with  experiment.  There  is  some  discrepancy  with  the  structured  grid  code  which  may  be 
on  account  of  the  fact  that  it  did  not  employ  inner  iterations.  Thus,  errors  arising  from 
factorization,  linearization  and  and  mismatch  of  operators  may  be  present  in  the  structured 
grid  solution. 

The  final  test  case  represents  an  exploratory  inviscid  computation  of  flow  over  a  multi¬ 
element  airfoil  system  undergoing  deployment,  as  an  example  of  bodies  in  relative  motion. 
The  geometry  represents  a  sectional  cut  of  the  wing  of  the  NASA  Langley  Transport  System 
Research  Vehicle  (TSRV),  and  was  obtained  by  direct  measurement  of  the  full  scale  aircraft 
(Boeing  737-100)  [46].  The  initial  position,  depicted  in  Figure  10,  corresponds  to  the  15° 
flap  setting,  while  the  final  condition  represents  the  40°  flap  setting,  which  is  the  approach 
configuration  for  this  high-lift  system.  The  initial  grid  about  the  5-element  airfoil  is  generated 
using  the  advancing  front  Delaunay  triangulation  method  [21]  at  the  15°  flap  setting  and  is 
displayed  in  Figure  10.  The  grid  has  26,191  vertices.  Figure  11  displays  the  Mach  contours 
for  steady  flow  at  this  setting  at  an  angle  of  attack  of  5°  and  a  freestream  Mach  number  of 
0.2. 

We  compute  the  time-accurate  flow  solution  to  full  deployment  using  a  five-grid  agglomer¬ 
ation  multigrid  procedure.  The  motion  is  prescribed  as  a  linear  variation  of  both  translation 
and  rotation  of  the  various  airfoil  elements,  and  is  computed  in  300  time  steps.  The  non- 
dimensional  time  for  deployment,  defined  as  t  =  Ij \[pool poo-,  is  taken  as  10.  Here  I  is  the 
chord  of  the  main  element  and  PooiPoo  are  the  freestream  pressure  and  density.  The  grid 
restructuring  involved  presmoothing,  spring  analogy  and  edge  swapping.  Figures  12  and  13 
depict  a  closeup  of  the  grid  in  the  flap  region  and  the  instantaneous  flow  solution  at  full  de¬ 
ployment.  Grid  quality,  while  not  as  good  as  with  the  initial  mesh,  is  acceptable  considering 
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the  large  scale  motion.  For  this  case,  rather  than  the  3-point  backward  approximation,  we 
have  chosen  to  use  the  implicit  midpoint  (trapezoidal)  rule: 

yn+l^rn+l  _  ynyrn  +  R{W^)  _  ^ 

When  the  connectivity  of  the  mesh  changes,  the  solution  changes  from  one  piecewise  linear 
representation  to  another  and  the  new  control  volumes  are  used  to  evolve  the  solution  to 
the  next  time  step.  Using  a  3-point  backward  rule  would  require  that  the  volumes  at  the 
previous  time  step  be  changed  as  well.  Figure  14  plots  the  total  lift  history  as  well  as  the  lift 
histories  of  the  elements  as  a  function  of  physical  time.  The  total  lift,  based  on  the  chord 
length  in  the  fully  nested  position,  increases  from  an  initial  value  of  2.53  to  a  final  value  of 
3.10. 

8  Conclusions 

An  efficient  implicit  time  integration  procedure  has  been  developed.  The  implicit  system  is 
solved  by  using  the  agglomeration  multigrid  procedure.  The  issue  of  the  mass  matrix  which 
arises  in  cell- vertex  methods  is  addressed.  It  is  shown  that  lumping  of  the  mass  matrix  may 
be  done  for  second-order  accurate  schemes  without  any  degradation  in  order  of  accuracy.  For 
higher  order  methods  based  on  point  values,  it  is  shown  that  the  lumping  of  the  mass  matrix 
degrades  the  order  of  accuracy  of  the  scheme.  Inviscid  and  viscous  calculations  have  been 
presented  for  flow  over  a  pitching  airfoil  and  an  impulsively  started  cylinder  and  the  results 
have  been  compared  with  experiments  and  other  computations.  A  mesh  point  movement 
strategy  has  been  proposed  and  tested.  This  has  been  used  to  compute  inviscid  flow  over  a 
multi-element  airfoil  system  undergoing  deployment. 
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Figure  9:  Histories  of  the  unsteady  residual  at  a  particular  time  step. 


Figure  10:  Lift  histories  during  the  third  cycle  of  motion  with  with  40,  20  and  10  steps  per 
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Figure  11:  Moment  histories  during  the  third  cycle  of  motion  with  with  40,  20  and  10  steps 
per  cycle. 
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Figure  12:  Moment  histories  on  the  coarse  0-mesh  without  swapping  of  edges. 


Figure  13:  Moment  histories  on  the  coarse  0-mesh  with  forced  swapping  of  edges 
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Figure  14:  Moment  histories  on  the  fine  0-mesh  with  forced  swapping  of  edges 


Figure  17:  Initial  grid  at  the  semi- deployed  position  (15°  flap  deflection). 


Figure  18:  Mach  contours  for  steady  flow  at  the  initial  position  (15°  flap  deflection). 
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Figure  19;  Closeup  view  of  the  grid  at  full  deployment  (40°  flap  deflection). 


Figure  20:  Mach  contours  of  the  instantaneous  solution  at  full  deployment  (40°  flap  defleo 


Appendix 

We  discuss  the  swapping  and  interpolation  procedures  in  three  dimensions.  Given  5  points 
in  three  space,  the  region  enclosed  by  the  convex  hull  can  be  tetrahedralized  in  one  of  three 
ways: 

1.  4  tetrahedra  sharing  a  common  interior  vertex. 

2.  2  tetrahedra  sharing  a  common  triangular  face. 

3.  3  tetrahedra  sharing  a  common  edge. 

Swapping  in  three  dimensions  consists  of  switching  between  the  second  and  the  third  choices 
to  improve  the  quality  of  the  tetrahedralization  [17].  Figure  25  illustrates  the  2-  and  3- 
tetrahedra  configurations. 
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(a)  2-tetrahedra  configuration  (b)  3-tetrahedra  configuration 

Figure  22:  Two  possible  tetrahedralizations  of  five  vertices. 

Conservative,  linearity-preserving  interpolation,  when  switching  from  the  2-tetrahedra 
configuration  to  the  3-tetrahedra  configuration,  is  similar  to  the  two-dimensional  situation, 
since  only  the  values  at  the  vertices  joined  by  the  new  edge  need  to  be  changed.  Define 

r2  =  ^  [(Wi  +  W2  +  W3  +  W4)Vi234  +  (VFi  -t-  W2  +  W3  +  W5)Vi245]  , 

T3  =  ^  [(Wl  +  W2  -f  W4  -t-  W5)Vi245  +  iW2  +W3  +  W4  +  W5)V2345  +  {W3  +  Wi  +  W4  +  W5)V3145]  • 

Changes  are  now  made  to  the  vertices  4  and  5  in  a  conservative  manner  by  distributing  the 
difference  T2  —  Tz  using  formulas  similar  to  Eqn.  (41),  where  the  W^ax  Wmin  taken 
as  the  maximum  and  minimum  over  the  five  points. 


Swapping  frona  the  3-  to  the  2-tetrahedra  configuration  is  different,  however.  Here  changes 
are  made  to  the  vertices  1,2  and  3  by  distributing  the  difference  Ts  —  T2.  Thus,  there  are 
two  degrees  of  freedom,  which  are  needed  to  enforce  the  condition  that  no  new  extrema  be 
created.  For  i  =  1,2,3  compute 


r,-  = 


T /nett;  Wmax—^Vi 
^  T3-T0  ’ 

T/nety  Wmi-n  -Wj 
T3~T2  ’ 


if  T3  r2  >  0 

if  T3  —  r2  <  0. 


and  define  ^ 

r  =  Min[ri,r2,r3,-]. 

Now  assume  r  =  rj.  Then 


=  Wiy/’"  +  r{T3  -  T2)/V^ 


Define  s  =  Min[r2,  ra,  0.5].  Changes  are  made  to  vertices  2  and  3  as  follows: 


ly/-  =  W2  +  ^(1  -  r) 

^neu,  =  W3  +  ^(1  -  r) 
6 


(5-r2)(a-r3) 

(0.5-r2)(0.5-r3) 


0.5  + 


(^-r2)(a-r3) 


(0.5-r2)(0.5-r3) 


0.5  + 


(3-0.5)(5-r3)  „  I 

(r2-0.5)(r2-r3) 


(^-0.5)(s-r2) 

(r3-0.5)(r3-r2) 


(1-s) 


(s— 0.5)(5— r3) 
(r2-0.5)(r2-r3) 


(1-5)  + 


($-0.5)(5-7-2) 
(r3— 0.5)(r3— r2)  . 


Formulas  for  the  remaining  cases,  when  r  =  5,^2,  ra,  can  be  derived  in  a  similar  man¬ 
ner. 
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